Astronomy & Astrophysics manuscript no. 5399pclt 


© ESO 2008 


February 5, 2008 





Kinematic frames and "active longitudes": does the Sun have a 

face? 

so"- 

O ! J. Pelt 1 , J.M. Brooke 2 ' 3 , M.J. Korpi 4 , and I. Tuominen 4 



Oh- 



1 Tartu Observatory, 61602 Toravere, Estonia 

2 Manchester Computing, University of Manchester, Oxford Road, Manchester, M13 9PL, UK 
rySs 3 School of Computer Science, University of Manchester, Oxford Road, Manchester, M13 9PL, UK 

■ 4 Observatory, PO Box 14, FI-00014 University of Helsinki, Finland 

m ; 

' ■ Received — ; accepted — 
y—i ■ 

j> . ABSTRACT 

in 

T^j- Context. It has recently been claimed that analysis of Greenwich sunspot data over 120 years reveals that sunspot activity 
£C) ■ clusters around two longitudes separated by 180° ("active longitudes") with clearly defined differential rotation during activity 
' cycles. In previous work we demonstrated that such effects can be observed in synthetic data without such features, as an 
artefact of the method of analysis. 

Aims. In the present work we extend this critical examination of methodology to the actual Greenwich sunspot data and also 
consider newly proposed methods of analysis claiming to confirm the original identification of active longitudes. 
^ C] Methods. We performed fits of different kinematic frames onto the actual sunspot data. Firstly, a cell-counting statistic was used 
O |' to analyse a comoving system of frames and show that such frames extract useful information from the data. Secondly, to check 
the claim of century-scale persistent active longitudes in a contramoving frame system, we made a comprehensive exploration 
of parameter space following the original methodology as closely as possible. 

Results. Our analysis revealed that values obtained for the parameters of differential rotation are not stable across different 
methods of analysis proposed to track persistent active longitudes. Also, despite a very thorough search in parameter space, 
• • ■ we were unable to reproduce results claiming to reveal the century-persistent active longitudes. Previous parameter space 
exploration has been restricted to frames whose latitudinal profile is opposite to solar surface differential rotation. Relaxing this 
restriction we found that the highest values of nonaxisymmetry occur for frames comoving with the solar surface flow. Further 
] ' analysis indicates that even these solutions are the result of purely statistical fluctuations. 
?d Conclusions. We can therefore say that strong and well substantiated evidence for an essential and century-scale persistent 
nonaxisymmetry in the sunspot distribution does not exist. 
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1. Introduction 

1.1. Nonaxisymmetry in late-type stars 

It is known that solar active regions may exist at the same 
position (such as recurrent sunspot groups) during several 
solar rotations, a new active region appearing at another 
position (e.g. Tuominen & Virtanen 119871 Pulkkinen & 
Tuominen 1 1 998f) . The sunspots reflect intense concentra- 
tions of magnetic field and thus occur only at certain dis- 
crete locations on the solar surface. For related discussion 
on the observed magnetogram recordings of the full-disk 
field see e.g. Stenflo IT9911 Consequently, different activ- 
ity indicators tend to be statistically correlated for a long 
time spans and a statistical analysis of the corresponding 
phase distributions must take this into account. 



The theory of the global magnetic field (mean-field dy- 
namo theory) has been developed over the last 50 years, 
starting from a linear theory (Parker 119551 Steenbeck 
& Krause 119691 for a detailed review see Krause and 
Radler 1980), followed by a nonlinear one, where the sta- 
bility of the solutions determines the symmetry proper- 
ties of the field (Krause & Meinel IT988I Brandenburg et 
al. ^989), and including a nonaxisymmetric solution as 
a possibility in the case of rapid rotation (e.g Radler et 
al. 119901 Moss, Tuominen & Brandenburg I199T|) . More 
recent models including hydrodynamics show that when 
the rotation velocity increases, the field geometry changes 
from antisymmetric with respect to the equatorial plane 
and axisymmmetric with respect to the longitude, to non- 
axisymmetric with respect to the longitude (e.g. Tuominen 
et al.HnH! Tuominen, Berdyugina & Korpi I2UU3|| . 
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Recently it has been claimed that long-term, i.e. persis- 
tent for at least 120 years, nonaxisymmetric structures can 
be detected in the solar sunspot data using new methods 
of analysis that take into account the possibility that such 
structures follow the differential surface rotation. This is 
an important claim and in the present paper we investi- 
gate it in detail because the existence of persistent "active 
longitudes" (nonaxisymmetric structures) in old slowly ro- 
tating late-type stars, such as the Sun, has important the- 
oretical implications. 

1.2. Nonaxisymmetry in solar data 

Berdyugina & Usoskin 200_3| (hence BU) analysed the lon- 
gitudinal distribution of regions of sunspot formation in 
the Greenwich sunspot data for a period of 120 years. By 
applying phase corrections, corresponding to the differen- 
tial rotation of regions of activity formation, they claimed 
to identify persistent (during 120 years) active longi- 
tude belts which are separated by 180° and which drift 
smoothly along longitudes. The essence of their method is 
to apply a longitude shift correction for each Carrington 
rotation to represent a different rotation velocity for the 
active longitudes as the sunspot activity wave moves from 
high to low latitudes in each solar cycle. As a response 
to the first paper (BU) Pelt, Tuominen & Brooke 2005 
(PTB) demonstrated that large part (or all) of the per- 
sistent smooth migration of the active longitudes can be 
the result of a particular method of data treatment in BU. 
Essentially it was shown that diagrams, similar to these 
used as a proof in BU, can be obtained even from random 
data. The criticisms of the paper were not directed at the 
observation that there can be active longitudes which per- 
sist over time, but towards the claim that such longitudes 
maintain consistent migration on a century time-scale. 

In recent papers Usoskin, Berdyugina & 
Poutanen l2"UfJ5l (UBP) and Berdyugina et al. |2"UUB1 
(BMSU) responded to this critique with the claim that 
they can confirm the original BU results using much sim- 
pler constructions which do not use the data processing 
methods criticized in PTB. In particular they consider 
that the new analysis confirms that the active longitudes 
maintain consistent smooth migration throughout the 
120 year period of the sunspot data analysed. The fact 
that this result was produced by two different methods of 
analysis of the data is considered to confirm the reality of 
the phenomenon. 

As will be discussed in length below, the migration 
models given in the first paper (BU) significantly differ 
from those obtained as a result of the analysis in UBP 
and represented in the data analysis part of BMSU. This 
is most readily seen in the cumulative frame kinematics: 
the BU frame makes about 28 extra rotations compared 
to the Carrington frame, contrasted with the UBP val- 
ues of about w 11 rotations for the Northern hemisphere 
and about ~ 14 for the Southern hemisphere. The par- 
ticular key feature here is that in BU the Northern and 



Southern frames keep coherently in phase whereas in UBP 
the Southern frame lags with respect to the Northern one. 
This indicates the need for a very thorough analysis of the 
methods of analysis used in constructing the differentially 
rotating frames considered to "freeze" the longitudes of 
these nonaxisymmetic features. 

This paper is structured as follows: in section [21 we 
describe the data and the background to our analysis. In 
section|21we demonstrate that kinematic frames which fol- 
low the rotation of sunpots exist and can clarify the basic 
properties of the longitudinal distribution of sunspots. In 
section 0] we present an examination of parameter space 
using the methods used in UBP but increasing the resolu- 
tion and the range of the parameter space search. In sec- 
tion El we discuss these results and introduce clarificatory 
statistical computations to check whether bimodal distri- 
butions could be statistical fluctuations. Finally sectional 
presents our conclusions. 

2. Data and methods of analysis 

2.1. Details of datasets analysed 

Our input data for statistical analysis was downloaded 
from the same Science at NASA web site 1 as in BU and 
UBP. The format description supplied with the data sets 
turned out to be incomplete and to unscramble the data 
we used the original document from NOAA Satellite and 
Information Service Site 2 . During the data checking we 
found some outliers in the data. The record from observa- 
tions of year 1980 which referred to longitude 408.9° was 
rejected and for couple of records where longitude was 
only slightly more than 360° we subtracted 360° from the 
given value (assuming circularity). The subset for years 
1878-1996 was then singled out to be compatible with pa- 
pers under scrutiny. In the terms of Carrington rotations 
we used rotations 325-1917 for both hemispheres. For each 
particular sunspot group or single spot we selected only 
the record of its first occurence as was done in original 
papers. It is important to note that some sunspot group 
numbers in data base are multiply used. To count them 
as separate entities we used the index value from the 21st 
column in the records (this was the piece of information 
missing in the first format description mentioned above). 
Finally we had two data sets: 1593 rotations and 18680 
records for the Northern hemisphere and 1593 rotations 
and 17966 records for the Southern hemisphere; the com- 
piled data sets are available on the web 3 . In the original 
papers "about 1600" rotations and "about 40000" records 
were used. 

To simplify the discussions below we will often renum- 
ber the original Carrington rotation numbers in our data 
set to start from 1 (hence 1, . . . , 1593). Furthermore, in all 
the plots and formulae we adopt the basic phase interval 

1 http:/ /science. msfc.nasa.gov/ssl/pad/solar/greenwch. htm 

2 ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/ 
SUNSPOT_REGIONS/GREENWICH/GROUPS/format.grp 

3 http:/ /www. aai.ee/~pelt/soft. htm 
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0° — 360°. As a result, some formulae can differ from those 
of the original papers. However, the original forms can eas- 
ily be recovered by using an appropriate phase shift. 

2.2. Construction of kinematic frames for analysis 

Two different physical scenarios can be considered when 
thinking of how to construct a framesystem. If one is to 
construct a frame to correct for the solar surface differ- 
ential rotation, it is natural to construct a frame that is 
moving with this flow. To investigate the situation where 
a rigidly rotating nonaxisymmetric structure anchored to 
the solar core would "illuminate" the solar activity belt 
influenced by differential rotation (the stroboscopic ex- 
planation discussed by BMSU), a more natural starting 
point would be to investigate frames which move against 
the solar surface differential rotation pattern. As both sit- 
uations could occur in the Sun, neither of them should be 
excluded from the analysis. 

Let us now divide the full set of sunspot records into a 
set of subintervals of length 27.2753 days, the Carrington 
rotation, as proposed originally by BU. For each particular 
rotation i,i — 1, . . . , N (where N is the total number of 
rotations in the set) the local phase shift fi^ is defined 
(in degrees per day). It can be constant or depend on a 
mean latitude of sunspot formation and free parameter 
B as in Eq. |SJ If the frame is to follow the differential 
rotation pattern then, over the full time span, the shifts 
accumulate according to the formula 

i 

Ai = A + TcJ2( n 3-^c), (1) 

i.e. the frame will lag if the sunspots rotate slower and will 
be advanced for more rapid rotation. Here Tc = 25.38, 
flc — 360° /Tq and A is a general phase shift to be com- 
puted later. Hereafter we refer to this frame system as 
the comoving frame. The other possibility is to choose a 
system of frames that moves in the opposite direction if 
compared to the surface rotation pattern, when the shifts 
are defined as 

i 

Ai = Ao + 2b5^(nc-nj). (2) 
i=i 

We will refer to this frame system as the contramoving 
frame. 

If the differential rotation is defined as Eq. HJ then we 
see that for the comoving frame, positive values of B lead 
to phases being pushed backwards if sunspots occur at 
higher latitudes where they rotate slower. If we allow for 
negative values of B, in the comoving frame they corre- 
spond to frames that are pushed forward with slower rota- 
tion. The opposite sign rule applies for the contramoving 
frames, positive Bs lead to phases being pushed forward 
with slow rotation and negative -Bs to a situation where 
the frame follows solar-type differential rotation. The basic 
difference of the two frame systems is illustrated in Fig0J 
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the comoving frame phase shift curve resulting from the 
best solution from the cell-counting statistics presented in 
Sect. |3 is plotted with a thin solid line and the original 
UBP Northern contramoving solution phase shift on top 
of it with a thick solid line. In the former solution the 
positive B-value causes the phases to lag behind during 
the early stages of the solar cycle (sunspots at higher lat- 
itudes) whereas the positive value in the latter solution 
results in phases being pushed forward during the early 
stages. 

If we do not restrict the parameter f^o and allow for 
both positive and negative values of the parameter B 
within each frame system, one can see that formally the 
two types of frames are equivalent, and both types of phys- 
ical situations of interest are detectable by both methods. 
In this study we utilise both types of frames with no re- 
striction on the sign of B. We note here that in UBP and 
BMSU the analysis was restricted to the contramoving 
formulation with positive B values only. It is a key point 
that, as we demonstrate below, although formally equiva- 
lent, they are not statistically equivalent. This is due both 
to the interaction of the frames with random fluctations 
in longitudinal distributions and also to the method of 
analysis proposed in UBP. 

To facilitate the comparison and transformation be- 
tween the framesy stems, the contramoving parameters 
translate to comoving frame by the following formulae 

j^comov j^contra 

Qcomov = _ Qcontra^ ^ 

Here we also note that special attention is needed when 
contramoving rotation laws are plotted together with 
comoving rotation laws (e.g. observational data from 
sunspots or helioseismology) . The contramoving param- 
eters must first be translated to comoving frame by the 
formulae given above. For instance, the best frame for 
the Northern hemisphere in UBP, obtained in the con- 
tramoving formulation with the parameters f2§°™ tra =14.33 
and B contra —3A0, translates to comoving parameters 
Clcomov^ u 05 and B comov = _ 3 40; j e slow ro t a tion at 

the equator which cyclically speeds up as a function of 
latitude. We have illustrated the possible error related to 
treating the comoving and contramoving frames as equals 
in our Fig. where we have plotted the best Northern 
UBP frame twice: the thin solid line represents the orig- 
inal illustration of UBP where the rotation law with the 
contramoving parameterization is plotted amongst comov- 
ing rotation laws without transformation. With the solid 
thick line we plot the transformed UBP rotation law in the 
comoving frame. Here it can be clearly seen that compar- 
isons of the UBP results with the solar surface differential 
rotation indicators are not meaningful. 

3. Cell-counting method 

To illustrate how the kinematic frames work in practice, 
in this section we construct and analyse one frame system 
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with a statistical method that is as simple as possible. 
We choose the comoving formulation, for which the phase 
shift is defined as Eq. ^ Next we subtract the flow model 
from actual longitudes 

Afei = Xki - A,; - mki x 360°, (5) 

and select integers mki so that the results will lay in the 
interval [0, 360°]. As we now work in the comoving formu- 
lation, after the subtraction we get a certain distribution 
of longitudes in a frame which is at rest with respect to 
the differential rotation pattern. 

Our comoving frame, as any other frame model in gen- 
eral, depends on a set of parameters (typically Oq, B and 
Ao). For each set of parameters we can compute corrected 
phases Xki (see Eq. |SJ and a value of a certain merit func- 
tion (criterion) which characterizes distribution of these 
phases (longitudes). For compatibility we will use merit 
functions whose parameter dependent minima are indica- 
tors of the desired types of phase distributions. 

To define our first merit fuction we divide all 
Carrington rotations into groups of ten rotations and all 
longitudes (phases) into 6° wide cells. We work with 1593 
separate rotations altogether and thus get a matrix of 
60 x 160 cells. For each set of free parameters (here Ao 
and f2o) we compute a merit function in the following 
simple way: Each of the corrected phases Xki belongs to a 
certain cell. Some of the cells remain empty because there 
are no sunspots in this particular rotation and interval of 
phases. As a criterion we use the ratio of occupied cells to 
the total number of phases. The rationale of this scheme 
is to produce a measure that is sensitive to the poten- 
tial existence of longitudes of enhanced spot formation. 
If the spots or spot groups form at persistent longitudes 
(in a corresponding frame), they tend to have their first 
observation in the cells which contain or are close to this 
persistent longitude. By fixing the cell length in number of 
rotations, we also fix the length of time interval for which 
we assume the activity to be persistent. The width of the 
cells fixes the amount of allowed phase dispersion during 
this number of rotations. 

3.1. Comoving frames with fixed angular velocity 

We now demonstrate that the formalism just described is 
sound and usable. We start from the simplest system of 
shifts f2j = flo, that is, we assign a constant shift f^o for 
each Carrington rotation. We let f^o vary in a wide range 
around the Carrington rotation velocity flc, and calculate 
the simple cell-counting merit function described above as 
function of fio; this is depicted in Fig. H] The dependence 
on Ao is "minimized out" . The curve is somewhat noisy, 
but still demonstrates clearly how the phases start to line 
up for Oo values larger than sa 14.1 deg/day and how the 
merit function value again increases from approximately 
14.6 degrees per day. The wide depression is too noisy to 
single out an optimal frame for all latitudes together. 

We now generalize our method to accomodate differen- 
tial rotation by preselecting latitudinal strips and repeat- 
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Fig. 1. The cell-counting merit function plotted against 
fio for the comoving frame with fixed angular velocity. 
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Fig. 2. Merit functions for the subsets of data in different 
latitude strips. The minima are now shifted and this is a 
result of differential rotation. 



ing the analysis on each of these subsets separately. We 
see from Fig. [21 that the minima for each subset are more 
clearly defined and show a shift representing the migra- 
tion between latitudes with different rotation velosities. 
The values of f2 for different strips can be used to build 
a latitudinally dependent rotation profile. The equatorial 
angular velocity is w 14.37, for middle latitudes w 14.18 
and for high latitudes ~ 13.40 degrees per day. The merit 
function curves are still noisy, but they can be improved 
by smoothing. Because of the different amount of activity 
indicators in different latitude strips the general level and 
scatter for the three curves are different. By constructing 
this merit function we can show that 

— The hypothesis of comoving frames shows the existence 
of a range of preferred frames even when applied to the 
full 120 year data set. 

— The scheme with phase corrections allows the com- 
putation of optimal parameters for different comoving 
frames. 

— The comoving frames in different latitude strips have 
stable rotational velocities revealed as principal merit 
function minima. 
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Fig. 3. Merit functions computed for the whole data set: 
B = (fixed angular velocity, thick line), B = 2.32 (vary- 
ing angular velocity, thin line). 

— The minima of the merit functions are located inside 
wide depressions, which shows that they are not minor 
fluctuations manifesting themselves as narrow peaks. 

— The kinematics of frames corresponds to kinematics of 
differential rotation - frames at higher latitudes tend 
to rotate slower than frames near equator. 

In this analysis we have used the longitudes of formation 
of all sunspots (not weighted by area) since we wish first 
to test the most basic hypothesis of persistence of activity 
formation. 

3.2. Comoving frames with changing angular velocity 

Moving on from the basic analysis above, we now attempt 
to detect an optimal single frame for all latitudes including 
weighting for sunspot area. Following UBP we will now 
include an additional term for each Carrington rotation: 

fi l = r> -Bsin 2 (V} l , (6) 

where (ipi)i is an area weighted average latitude for i-th 
Carrington rotation and parameter B measures the am- 
plitude of the latitude dependent part. 

To get the optimal triple of frame parameters we need 
to compute merit functions for a large grid of trials. Using 
a two step procedure, decribed in detail below, we found 
a global minimum for the parameter grid (flo x AqX B) = 
[13.5 - 15.0; 0.001] x [0.0 - 6.0; 0.02] x [1 - 5; l] 4 to be 
located at n Q w 14.416 and B « 2.32. The range of A 
investigated results from the fact that there are only six 
different starting positions for 6°-wide cells. The resulting 
merit functions are depicted in Fig. [3| and the correspond- 
ing shifts for the comoving frame in Fig. 0] 

It is important to stress that the strongest minimum 
for the frames with changing angular velocity is deeper 

4 Here and below we use a systematic notation for compu- 
tation grids. Inside the square brackets we give minimum and 
maximum values for the parameter in question followed by the 
stepsize of computation. 




400 500 600 700 800 900 1000 1100 1200 1300 
Rotation 

Fig. 4. Typical shift curves for comoving (thin line) and 
contramoving (thick line) frames. 

compared with that of the fixed angular velocity frames. 
Consequently it indeed models to some extent real features 
in the distribution of sunspot longitudes. 

This analysis demonstrates that the formalism of co- 
moving frames can reproduce previous results which show 
that sunspots form at longitudes which can persist for time 
period much longer than a Carrington rotation. 

4. The UBP method 

The basic claim of BU, UBP and BMSU is that two pre- 
ferred longitudes separated by 180° persisting over the last 
120 years can be seen in the distribution of sunspot lon- 
gitudes. These "active longitudes" migrate in any frame 
with fixed angular velocity, but form coherent and persis- 
tent longitudinal structures in a certain frame with chang- 
ing angular velocity. In the previous section we have shown 
that frames with changing angular velocity can be used to 
reveal short-term (of the order of 10 rotations) persistent 
structures. 

The first difference between our simple cell-counting 
example method and the method of UBP comes from the 
fact that they have chosen to use phase corrections of the 
contramoving type, given by EqEl The second difference 
is the usage of sunspot areas. Our simple merit function 
treats only longitudes whereas the areas are accounted for 
only when computing the mean latitudes. In UBP the ar- 
eas are first normalized and then used as weights in the 
merit function. This amplifies the effect of the spots oc- 
curring during the low activity interval of the solar cycle. 
The third difference is the usage of different merit func- 
tions. Since in UBP the authors seek a particular phase 
distribution they use a particularly crafted bimodal merit 
function (see Fig. Eland Eq.0, whereas ours seeks to min- 
imize phase drift between cells. 

Since the claims of UBP go well beyond previous anal- 
ysis of persistence of nonaxisymmetric features and lead 
to results which are difficult to explain (e.g. disconnection 
of Northern and Southern hemispheres) we present a very 
detailed and careful analysis of the full parameter space 
below. The only principal difference between our analysis 
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Fig. 5. Comparison of differential rotation curves ob- 
tained from helioseismology (SOHO/MDI, Schou et 
al. H.998), for sunspots of different ages (Pulkkinen & 
Tuominen 11998(1 . and for the different frame formulations 
investigated in this paper. 



in this section and the UBP investigation is the extension 
of the method to negative values of B, i.e. we also allow 
for the possibility to correct for solar-type differential ro- 
tation. Otherwise we have followed the original implemen- 
tation as carefully as possible. The important refinements 
are given in detail below. Possible reasons why our results 
differ from those of UBP are given in Sect 03 

4.1. Search for the optimal frame 

For every combination of parameters f^o, B and Ao we can 
build a corresponding frame. The shifted phases arc then 
evaluated with the UBP method using the merit function 



ik 



(7) 



where = min(min(Aifc, 360° — Xik), \Kk ~ 180°|) mea- 
sures the distance between the corrected phases and the 
nearest centre (0° or 180°) 5 . 

In the actual computations described below we used a 
phase-shifted version of the merit function, so that the ex- 
pected maxima in the longitude distribution were shifted 
to 90° and 270°. This is to ensure uniformity of all the 
plots. 

Thus our absolute phase values differ from UBP, how- 
ever estimates of bimodality and departure from axisym- 
metry are clearly unaffected by this shift. In the plots 
below we depict the merit function values without nor- 
malization, this does not affect the interpretation of the 
plots. In the tablen we present the true (normalized) val- 
ues of £. 



5 In the printed version of UBP this function is given incor- 
rectly. This mistake was also confirmed by I. Usoskin, one of 
the authors of UBP I5U051 
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Fig. 6. The merit function depends on the distance of a 
particular phase from the nearest centre (90° or 270°). It 
is expected that the distribution sought for will have two 
maxima, situated 180° from each other. 
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Fig. 7. The merit function dependence on Slo is computed 
with two different time steps: 0.01 and 0.001. 



The dependence of Cik-s on the corrected phases is de- 
picted in Fig. [S] Values of Ay. in Eq. [7| are defined as 



A-ki — Ski I ^ ] Sji, 



(8) 



where S is the observed spot area corrected for the pro- 
jection effect, and the sum is taken over all the accounted 
for spots in the given Carrington rotation. 

For each particular triple (Ao, B) we can compute 
the corresponding shifts and then evaluate the obtained 
phase (longitude) distributions using the merit function 
£. Before the actual computations for a full parameter 
space it is important to estimate the required step lengths 
along each parameter. 

We start from the parameter Uq. If we look how the 
shifted phases depend on changes in this parameter, we see 
that a small incremental change in the parameter, Af^ is 
multiplied by the constant Tc = 25.38 and then, depend- 
ing on the particular rotation, by the number of accumu- 
lated rotations. For the last phase the number of rotations 
is approximately 1600. 
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Fig. 8. The merit function dependence on B is computed 
with two different time steps: 0.02 and 0.001. 

To recover full details of the merit function we fix a 
certain maximum allowed phase shift for one step change 
in the input parameter Slo- Because there are two increas- 
ing and two decreasing parts in the distance function (see 
Fig- El 7 we choose one eighth of the full phase range (45°) 
as the maximum allowed phase shift. We obtain 

Afi Q = 45/(T c * 1600) « 0.001. (9) 

The suitability of this choice is clearly illustrated in Fig. [3 
where we plot the merit function running along f^o with 
two different steplengths: 0.01 (thick curve, the stated pre- 
cision in UBP) and 0.001 (thin curve). Some important 
features are certainly lost for the undersampled case. (In 
this plot and in all the following or B dependent plots, 
the dependency on Ao has been "minimized out" using the 
steplength 5°.) For local refinements and for short param- 
eter intervals we have used even smaller values of Af2o- 

The precise analysis for the B parameter is complicated 
because of the random component in the mean latitudes. 
Some simple conclusions can be made by using a trial and 
error analysis. As seen from Fig. [S] the choice of the step 
length 0.02 (thick curve, stated precision of B in UBP) 
is reasonably good, compared with the steplength 0.001 
(thin curve). Nevertheless, some minima can be slightly 
misplaced due to the undersampling. As a compromise, a 
proper choice for the B stepsize could be 0.002. 

The dependence of the merit function on the parameter 
Ao is obviously smooth and we can use a coarser step in 
degrees, say 5° to 10°. Because the distance function is 
periodic with a period 180° we can restrict our search to 
the subinterval 0° - 180°. 

We can now compute the total number of trial triples 
we need to work with. Say we want to seek for a global 
merit function minimum in the parameter space [10 — 
20] x [0 — 6] x [0 - 180] (according to one of the authors - 
I. Usoskin l|2005[) . this was exactly the chosen parameter 
range in their initial search). With the chosen steplengths, 
we need to evaluate at least (10 * 1000) * (6 * 500) * 
18 = 540000000 triples. This is far beyond our computa- 
tional capabilities. Computations with significantly longer 
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Fig. 9. The merit functions for three different values of B. 
It is clearly seen that all three are quite similar and only 
slightly shifted along f^o- On a two-dimensional diagram 
this would look like a row of slanted stripes. 

steplengths for the parameters would leave us uncertain 
of recovering the actual global minima. However, we are 
able to reduce the search space by observing that B and 
flo are not completely independent parameters. Observing 
the similar shapes of curves for different values of B in 
Fig. El and noting the structure of Eqs. and El we see 
that the cumulative shifts depend monotonically on the 
positive-definite sums 

i 

s, = J2 sin2 Mi- ( 10 ) 

i=i 

This component effectively consists of two parts: a wave 
which models the differential rotation during the solar cy- 
cle and a linear part which depends monotonically on the 
sum of previous waves. In the final cumulative sums the 
fixed part of this sum is combined with a similar linear 
component which comes from the parameter Slo- Because 
these contributions are included in the cumulative sums 
with opposite sign, we can compensate for the increase 
in the parameter B with the increase in the parameter 
Qq. Their total contribution, which is the difference, then 
stays fixed. The only change is a small redistribution in 
phases which comes from the cyclical part. 

We can now reparameterize our analysis. First we de- 
fine a linear part of the curve Sj, j = 1,...,N. Then we 
subtract this linear part from the curve and add it to the 
fio dependent part. Finally we will have new parameter 
Q' which will depend on original Qq and B: 

n' = n Q -A c x b, (ii) 

where Ac = {Sn — Si)/(N — 1) is the linear slope of 
the cumulative sum of Sj-s. The new parameter tt' Q will 
now describe an overall constant rotation incrementing 
in phase and the parameter B measures the contribution 
from the differential rotation. The parameter Sl now de- 
scribes the mean comoving framespeed and the parameter 
B small phase shifts due to the differential rotation. The 
effect of reparameterization is clearly seen in Fig. ^1 the 
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Fig. 10. For the reparameterized case the dependence on 
B is much smoother (compare with Fig.[S|). 



Fig. 11. The merit function computed with B = for the 
Northern hemisphere. 



dependence on the parameter B is much smoother and it 
is possible to decrease the required number of parameter 
triples. We can now formulate our optimization strategy. 
First we perform a coarse search using the reparameter- 
ized triples. The proper steplengths for this case are 0.001 
for fig, 0.1 for B and 5° for Ao. The standard tto can then 
be computed from Q' . In the vicinity of the detected min- 
ima we can perform high resolution refinements using the 
normal triples of parameters with steplengths 0.0001 for 
fl , 0.002 for B and 1° for A . 

To characterize the search results and estimate their 
significance we adopt the measure of nonaxisymmetry de- 
fined by UBP. For each distribution of longitudes we form 
two subsums of corresponding spot areas 



Ni = Y, A ^, if |A H - 90°| < 45° or |A fel - 270°| < 45°, 

k.i 

N 2 = A ki , otherwise, (12) 



k.i 



where the summation is taken over all the spots in all 
the Carrington rotations. The nonaxisymmetry T is then 
defined in UBP as 



r = 



Ni -N 2 
N!+N 2 ' 



(13) 



Non-uniform (in longitude) distributions give higher val- 
ues of r. 



4.2. Search results 

We restricted our search space for fio to the interval 
[13.5 — 15.0] since this includes all the physically plau- 
sible values for this parameter. We present our results in 
three parts. First we give results for fixed angular velocity 
frames, then for frames with particular values of B and 
finally for the full parameter space investigated. This al- 
lows us to illustrate some important aspects of the merit 
function "surface" . 



4.2.1. Fixed angular velocity frames 

From the previous discussion we know that it is quite easy 
to find particular comoving frames for separate latitude 
strips and even mean frames for the full data. Therefore 
we first evaluate frames with B = 0. From UBP we learn 
that "...the hypothesis of rotation of the active longitudes 
with a fixed rate (giving T = 0.02 — 0.03; see Figs. 2a, 
3a, and 7a,b) cannot be distinguished from the null hy- 
pothesis of the axisymmetric sunspot distribution" . Our 
results of the subset search with B = are depicted in 
Fig. ITTl and in tablc^ Comparing Fig.^^with Fig.^ we 
observe a key difference: instead of the wide depression 
we see a number of sharp minima. The strongest peaks 
are located at « 14.05 for North and « 13.85 for South. 
The nonaxisymmetry measure is T — 0.064 (North) and 
T = 0.073 (South), instead of 0.02 - 0.03 in the original 
paper. We are dealing here (and also below) with the fluc- 
tuations themselves, therefore error bars (e.g. about some 
mean value) are not relevant. 

4.2.2. Frames with fixed values of B 

To make straightforward comparisons with the UBP re- 
sults we sought merit function minima of slices with fixed 
B values. For the Northern hemisphere we tried to recover 
the original result in a slice with the best UBP value of 
B = 3.40; for the Southern hemisphere our slice was se- 
lected with the best UBP value of B = 3.39. The results 
are given in tabled 

In addition to the computations with positive values of 
B we computed "spectra" also around their negative coun- 
terparts. It is clearly seen from the corresponding tables 
that our results differ significantly from these presented in 
UBP. Fig. E| is of particular significance when compared 
with Fig. |21 where we can see that the best comoving frame 
shows itself as a very clear minimum in the middle whose 
width is much greater than that of the fluctuations. In 
Fig. 1121 we observe only a bunch of narrow fluctuations. 
With the merit function used in section|21 when we change 
the parameter values, then in the vicinity of the main min- 
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Table 1. Summary of the search results with the UBP method. The analysis is performed in the contramoving 
formulation: the obtained frame parameters are given in the columns 5 and 6, and transformed into the comoving 
formulation in columns 3 and 4. The corresponding best values from the original UBP paper, if available, are given in 
columns 9 and 10. 
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Fig. 12. The merit function computed with B = 3.40 for 
the Northern hemisphere. This figure should be compared 
with Fig. 



imum, phases migrate slowly from one cell to another and 
the general picture changes slowly. With the UBP merit 
function (crafted for bimodal waves), we observe a highly 
fluctuating curve, with no observable minimum of width 
greater than the fluctuations. 



Fig. 13. Strongest fluctuation for Northern hemisphere. 



= -5.10 Qo - 14.120 



4.2.3. Global search 

In the global search we allow B also to vary freely. The 
results of the global computations are summarized in ta- 
ble n and in Fig. El The distribution for the Northern 
hemisphere does not show the clear bimodal structure pre- 
sented in UBP, but its nonaxisymmetry is quite high - 
r = 0.096. Our results with the contramoving formula- 
tion and UBP merit function show that there are no min- 
ima with width greater than the fluctuations which would 
be expected from causal dependence on the parameters. 
Given this, we observe that the strongest minima tend 
to occur for the negative values of B, which translate to 
solutions with the comoving frame formulation. 

Possible reasons for the difference between our results 
and the results presented in BU and UBP can be: 



50 100 150 200 250 300 350 

Longitude (deg) 

Fig. 14. Longitudinal distribution obtained with best pa- 
rameter fit for a merit function with four maxima. The 
parameter space used in the search was the same as in 

Fig.ua 



— The data sets used in UBP and here can be different. 
The authors of UBP may have used certain pre- or 
postselection principles applied to the input data (not 
described in the paper). 

— The curves of the mean latitudes (and therefore also 
the shift frames) in UBP were smoothed in a particular 
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way (see Fig. 1 in UBP) or were interpolated in a par- 
ticular way (compare Fig. 5 in BMSU to our Fig. I15|) . 
However, we did not find any comments about these 
procedures in either of the papers. 

Below we try to look at some other possibilities. To check 
if the bimodal merit function gives clearly superior results 
for departures from axisymmetry, we constructed a merit 
function with four maxima over the full range of longitude 
and performed a full parameter space search as above. 
This represents not two migrating density waves but four, 
two for each side of the Sun. The fluctuation which gave 
the global minimum for this criterion is depicted in Fig. 1141 
(Northern hemisphere) . The distribution is certainly non- 
uniform and four expected waves are clearly seen. 

5. Discussion 

The results obtained above and those of BU, UBP and 
BMSU (data analysis part) show discrepancies which re- 
quire some explanation. Below we present how we tried 
to control our procedures to discover potential reasons for 
such major differences. We also include comparative anal- 
ysis of the results of BU and UBP since a major contention 
of PTB is that the attempt to detect persistent features 
of the sunspot data is very sensitive to the methodology 
used. 

5.1. Kinematics of frames 

The most significant problem with the persistent longitude 
strips is connected to their cumulative kinematics. 

The authors of BU found that in the Carrington refer- 
ence frame the observed migration had a phase lag about 
2.5 solar rotations per sunspot cycle, in total about 28 
rotations during 120 years. 

In the follow-up study (PTB) we demonstrated that 
a major part of this effect can be due to hidden assump- 
tions in the methodology of data processing. In the next 
paper (UBP) the authors attempt to confirm their origi- 
nal conclusions by constructing a new set of contramoving 
frames using global optimization techniques. The selected 
dynamic frames, however, differ significantly between BU 
and UBP & BMSU: in the former case the pattern makes 
28 extra rotations during 120 years, whereas in the lat- 
ter, the Northern pattern makes w 11 and Southern « 14 
rotations. The reader can easily verify this by compar- 
ing Figs. 3-4 from BU with Fig. 1 in UBP and especially 
with Fig. 1 in BMSU, where the different frames and their 
kinematics are depicted graphically. Consequently we con- 
sider that the two latest papers (UBP and BMSU) fail to 
confirm the results obtained in the first paper (BU), in 
particular they do not refute the conclusions of PTB. 

The second kinematical problem is connected to the 
different parameterization schemes. In the UBP investi- 
gation, a contramoving frame system is used, in addition 
restricted to positive values of the parameter B. As dis- 
cussed in length in Sect. 12.21 to construct a frame sys- 



tem comoving with the solar surface differential rotation 
pattern, negative B values should have been allowed for. 
Now the solutions are restricted to trace a situation where 
movement in the opposite direction than the differential 
rotation pattern occurs. Since we have demonstrated in 
section l2~2l that these frames are not equivalent methods 
of treating the observational data, the restriction of pos- 
itivity of B should be dropped, as in this study. As a 
consequence, with the UBP method, we find the highest 
nonaxisymmetry measures for frames with negative Bs, 
i.e. moving with the solar surface differential rotation pat- 
tern. 

The third kinematical problem comes from the fact 
that the numbers of extra rotations (compared with 
Carrington frame) are in the UBP model significantly dif- 
ferent for the Northern and the Southern hemispheres: 
w 11 and w 14 correspondingly (cf. Fig. 1 in BMSU). 
This raises serious difficulties for the dynamo theory of 
the large scale magnetic field. 

5.2. Evidence of flip-flops in the dominant longitude 

Another important result of BU is the alternation of dom- 
inance between the two longitudinal maxima with the ap- 
proximate periods of 3.8 (North) and 3.65 (South) years 
(periodicity in flip-flops). The flip-flop phenomenon was 
originally found in an extremely active late-type star FK 
Com (see Jetsu, Pelt & Tuominen ll993fl . 

The particular values for the mean flip-flop periods for 
the Sun (in BU) were obtained from an analysis of the 
power spectra. Firstly the evidence presented does not es- 
tablish that such periodicity exists. From the power spec- 
trum plot of Fig. 8-9 (BU) it is clear that nearly all the 
spectral power fluctuations are above the line of 95% con- 
fidence level. If there were clear and dominant periodicity 
one would not expect this. 

In UBP it is claimed that the flip-flop events and 
their periods obtained in BU are valid also for the newly 
constructed dynamic frames and the resulting coherent 
nonaxisymmctric structures. The flip-flop event frequency, 
however, was computed using a frame which rotates nearly 
three times more rapidly than the dynamic frames found 
in UBP. We note that a phase shift of ps 180° which causes 
the flip-flop event position to swap longitude can be pro- 
duced by 0.004 degrees per day change in the frame model 
parameter f^o (quoted precision in UBP for this parameter 
is 0.01). In the new construction, in addition, the shifts are 
cumulative sums of functions from a significantly fluctu- 
ating mean latitude curve (see Fig. 115(1 . Thus evidence for 
the persistence of this phenomenon across the two meth- 
ods of determining rotating frames is not sufficiently es- 
tablished. 

Consequently it is not reasonable to combine the re- 
sults obtained from one frame system with the results of 
a totally different system. 
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Fig. 15. The mean latitude of sunspots (weighted by 
sunspot area) against time. The magnitude of the fluc- 
tations against the general trend can be seen. 



5.3. Nonaxisym metric distributions in fixed velocity 
frames 

It is a major point of the UBP analysis that measures of 
nonaxisymmetry are low (r = 0.02 — 0.03) for frames with 
fixed rotation (see Fig. 7 in UBP). In our computations 
we obtained a relatively high T for fixed rotation frames 
(see Table nj, of the same order as the best values for 
differentially rotating frames. Since we have demonstrated 
that the UBP merit function shows rapid fluctuation with 
respect to changes in the parameter values, it is also very 
possible that results are sensitive to preprocessing of data. 
Accordingly we have published on the web the exact data 
that we used (c.f. section[^J). We think that in looking for 
evidence of persistent longitudes it is necessary to examine 
parameter space very carefully and be completely explicit 
about all data processing steps used, with publication of 
data as an aid to other researchers. 



5.4. The effects of folding phase shifts 

In the process of seeking other possible reasons for dif- 
ferences between our computations and those of UBP we 
folded (by 360°) the shift curves with different parame- 
ter values. In Fig. ^] the shift curve for the UBP con- 
tramoving frame is depicted. The key features in this plot 
are the "S"-shaped fragments. For different frame param- 
eters these fragments line up differently. The correspond- 
ing marginal distributions in Fig. 1171 show that the shifts 
themselves are not distributed evenly. In Fig. ^] we have 
depicted how the merit function changes for a fixed value 
of B — 3.40 and changing value of Slo- The plot demon- 
strates that the merit function has a specific wide depres- 
sion; in the bottom of this depression around f^o values 
14.20 — 14.40 there are deeper minima. Fig. ^] (inlet) re- 
veals that the deepest minimum occurs at fl = 14.316. 
The second deepest is at 57 = 14.330. Comparison of 
Fig. U3 and Fig. d shows that the value fi = 14.330 is 
the strongest peak whose shift system has a characteris- 
tic double mode distribution. Just to remind the reader, 



Fig. 16. Particular shifts for a contramoving frame with 
n = 14.330 and B = 3.40 folded with module 360°. 
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Fig. 17. The marginal distribution (along longitudes) of 
the frame depicted in Fig. 1161 
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Fig. 18. The merit function computed with longitudes set 
to zero. There is a strong depression around fig = 14.20 — 
14.40. Enlargement of the depression part is given as inlet. 

we did not use longitudes at all to compute the last three 
plots. We can see how sensitive the form of the distribu- 
tion is to parameter values giving merit function minima 
of comparable strength. 

We draw also attention to the particular value of 
Qo = 14.33 deg/day, which occurs both as a solution in 
UBP and as fluctuations in the input data (shift curve) for 
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Fig. 19. The marginal distribution for the frame with 
O = 14.316 and B = 3.40. (compare with Fig. lT7jl . 



the Northern and Southern hemispheres: is there anything 
peculiar behind this coincidence? The unexpected solu- 
tion can be found with the help of the previous analysis. 
To compute minima for the large parameter grid we used 
certain reparameterization where we divided the frame 
shift curves into two parts: the linearly rising part and 
the cyclical part. From the linear part we computed then 
mean rotation velocity for the frame. For the particular 
values B = 3.40 and fi = 14.330 this gives % = 14.094. 
The charming magic of this value reveals itself when we 
compute first how much it differs from the Carrington 
siderial value flc — = 0.090. From the angular velocity 
we can compute the corresponding period which occurs 
to be 10.96 years - closely coinciding with the mean so- 
lar sunspot cycle length. The obtained result can also be 
read out from Fig.l in BMSU - there are approximately 
11 cycles and the range of the shifts is also approximately 
11 full rotations. 

5.5. Selection of a correct null hypothesis 

We now consider the selection of a null hypothesis to check 
the significance of evidence for persistence in the longi- 
tudes of sunspot activity. In UBP the significance of the 
obtained results was evaluated using a reference (null hy- 
pothesis) distribution which was computed using Monte- 
Carlo type methods. For each Carrington rotation they 
randomly permuted all the sunspots, i.e., a new set of ran- 
dom Carrington longitudes were ascribed to each actually 
observed sunspot while keeping its area. Then the value 
of r was calculated as described above. They computed 
the nonaxisymmetry T for 5000 sets of such random-phase 
sunspot occurrences to get a reference distribution. This is 
a valid procedure for the case when we need to test certain 
particular longitude distributions against totally random 
distributions. 

However, this analysis does not take into account the 
fact that we select the best fit solution among number 
of candidate solutions. Consequently we need to compute 
the "false alarm probabilities" (see for detailed explana- 
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data. 
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Fig. 21. Distribution of records for short-term correlated 
data. 

tion Scargle 1982). Thus in the Monte-Carlo framework 
we must compute distributions for "record values" . The 
full scale, full resolution Monte-Carlo computation for full 
Qo,B parameter space is certainly out of our computa- 
tional capabilities. Fortunately, the general picture be- 
comes evident from a subset of calculations. First we re- 
stricted our computations to parameter space slices with 
fixed parameter B = 3.40 and a grid with moderate res- 
olution for ilo - [13.5,15.0,0.01] (Northern hemisphere). 
Then we performed for 1000 randomized (as in UBP) data 
sets search for strongest peak in Qq spectrum. For the 
strongest peak we computed the value for T. The corre- 
sponding value distribution is depicted in Fig. [3UJ It is 
evident that, even for fully decorrelated data, the record 
values can be quite high. This can be explained by the 
fact that all data points are not equal - they are weighted 
by area. The areas differ very significantly (from to 
1894 millionths of total area) and the number of effec- 
tive degrees of freedom to be distributed along longitudes 
is significantly less than total number of activity elements 
(18680). 

It is known that in the Sun we can observe quite long- 
lived activity complexes, for instance hot spot systems of 
solar flares with lifetimes of the order of the solar cycle 
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length (BaiUHHlEnnni)- To model this kind of situation 
we performed an additional test. Instead of full random- 
ization of the longitudes we divided the full observational 
data into 10 rotation length sections. For each of the 1000 
Monte-Carlo runs we added random shifts from the in- 
terval [0, 360] to every section. In this way local correla- 
tions remained untouched, but long range correlations (if 
present) were certainly washed out. The record distribu- 
tion for this experiment is depicted in Fig. [2] and it is 
to be compared with results in table ^ We see that we 
can safely treat global minima as the random fluctuations 
due to the local correlations and the uneven distribution 
of sunspot areas. Our analysis is conservative - we under- 
sampled along fl parameter, we used a fixed value of B 
and we did not seek minima in T but in the UBP merit 
function. 

In addition to the above mentioned correlations there 
is another source of flexibility in the UBP method which 
results in additional fluctuations. This is due to the 
reweighting procedure chosen. In the input data set there 
are some rotations which contain nearly 70 new activity el- 
ements and all these get minuscule weights in the method. 
The weights for the elements which happen during the ro- 
tations between the activity cycles, on the other hand, 
are very powerfully amplified because these rotations are 
strongly underpopulated (rotations with only a single new 
event occur reasonably often). As a result, many kinds of 
distribution may form (see for instance Fig. I14f). Instead 
of amplifying the effect of persistent activity waves the 
method amplifies random fluctuations. 



5.6. Inherent constraints of the data set. 

Given the difficulties we have described above we examine 
the constraints caused by the limited time span of the data 
set, relative to the solar cycle. The normal comoving flow 
of the activity indicators can be appropriately measured 
by such methods. But for activity waves which travel in 
the opposite direction to solar rotation the evidence can 
only be quite a weak difference between a low level real 
physical effect and the fluctuations which result from the 
statistical nature of the data at hand. 

As shown in PTB it is possible to combine a small 
number of long stretches of data (using appropriate phase 
corrections) so that their marginal distribution is almost 
entirely an artefact of the method of analysis. 

The inherent (for such an analysis) phase ambiguity is 
a strong constraint which can be overcome only by using 
much longer (in time) data sets. The current time extent 
of the sunspot database covers only a small number of 
correlation length size subparts and this results in a high 
level of random fluctuations. We do not deal here with 
1600 independent rotations but with a significantly lower 
amount of activity complexes. 



6. Conclusions 

In the recent papers BU, UBP and BMSU the authors 
claim that, hidden in the sunspot distribution, there are 
two persistent active longitudes which migrate according 
to a certain differential rotation law. After presenting a 
critique of the BU results in the earlier paper (PTB), here 
we have tried to check the results obtained in UBP and 
the data analysis part of BMSU. Our results can be sum- 
marized as follows: 

— The contramoving frames constructed by BU differ sig- 
nificantly from the frames found using mathematical 
optimization techniques in UBP and BMSU. In the 
first case the frame does approximately 28 extra ro- 
tations (compared with the Carrington frame) and in 
the new papers the number of extra rotations is around 
10-11 for the Northern hemisphere and 13-14 for the 
Southern hemisphere. Consequently the claim in UBP 
that the new analysis confirms the results obtained in 
the previous study is not well founded. 

— In the two latest papers authors re-iterate the claim 
presented in BU of the presence of the flip-flop effect 
with a mean period of 3.7 years. This result was ob- 
tained from (and now presumably incorrect) analysis 
in the first paper (BU). 

— Our analysis demonstrated that even for fixed angu- 
lar velocity contramoving frames it is quite easy to get 
particular longitude distributions with relatively high 
values of nonaxisymmetry. This result is quite robust 
and does not depend on the particular implementation 
of the UBP method (smoothing of the shift curve, sam- 
pling rates in the parameter space etc.). Consequently 
all the claims in UBP about T = 0.02 - 0.03 level fluc- 
tuations in the case of B = arc incorrect. The pa- 
rameter space is full of parameter combinations which 
can give a rather high values for nonaxisymmetry. 

— Our re-implemantion of the UBP method gave funda- 
mentally different results from those of the original au- 
thors. We can only conjecture about the possible rea- 
sons for the discrepancies. The pre- or postselection of 
data points, parameter undersampling in search pro- 
cedures, incorrect phase reduction etc. - these all can 
be considered as possible sources of differences. 

— From inspection of the strongly fluctuating merit func- 
tion curves and additional Monte-Carlo type computa- 
tions for parameter space slices, we conclude that the 
frames obtained from the optimization procedure are 
pure fluctuations. 

— We found that the parameter values obtained as re- 
sults in UBP can be seen already in the distributions 
connected to the frame shift curve itself (which does 
not depend on longitudes of sunspots at all). 

— Our analysis does not indicate that short-lived features 
concentrated in certain longitude are totally missing 
in the sunspot distribution, even though evidence for 
century scale persistent active longitudes can not be 
found. The simple cell-counting method investigated 
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here merely as a demonstration tool may well be used 
to estimate real coherence times of these structures. 
— Whichever method used, the preferred solutions oc- 
curred at positive/negative values of the differential 
rotation parameter B in the comoving/contramoving 
frame systems, respectively, i.e. the frames had a ten- 
dency to follow a solar-type differential rotation pat- 
tern. 

From this we conclude that the results obtained in BU, 
UBP and BMSU (data analysis part) are inconsistent, and 
evidence for well established persistent activity migration 
is still lacking. In other words we are inclined to think that 
the enigmatic face of the Sun is still hidden and that the 
patterns which the authors of the three papers saw have 
a human origin. 
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